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Abstract 

We discuss the transient and steady state fluctuation relation for a mechanical 
system in contact with two deterministic thermostats at different temperatures. 
The system is a modified Lorentz gas in which the fixed scatterers exchange 
energy with the gas of particles, and the thermostats are modelled by two Nose- 
Hoover thermostats applied at the boundaries of the system. The transient fluc- 
tuation relation, which holds only for a precise choice of the initial ensemble, is 
verified at all times, as expected. Times longer than the mesoscopic scale, needed 
for local equilibrium to be settled, are required if a different initial ensemble 
is considered. This shows how the transient fluctuation relation asymptotically 
leads to the steady state relation when, as explicitly checked in our systems, the 
condition found in [D.J. Searles, et al, J. Stat. Phys. 128, 1337 (2007)], for the 
validity of the steady state fluctuation relation, is verified. For the steady state 
fluctuations of the phase space contraction rate A and of the dissipation function 
$7, a similar relaxation regime at shorter averaging times is found. The quantity 
Q, satisfies with good accuracy the fluctuation relation for times larger than the 
mesoscopic time scale; the quantity A appears to begin a monotonic convergence 
after such times. This is consistent with the fact that Q, and A differ by a total 
time derivative, and that the tails of the probability distribution function of A are 
Gaussian. 

1 Introduction 



In the last decades it has been recognized the central role that the fluctuations 
around stationary states play for a microscopic theory of nonequilibrium phenom- 
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ena [1]. Furthermore, biological sciences and nano-technology are ever more 
interested in nonequilibrium fluctuations. 

The 1993 paper by Evans, Cohen and Morriss [2], on the fluctuations of the en- 
tropy production rate of a deterministic particle system, modeling a shearing fluid, 
provided a unifying framework for a variety of nonequilibrium phenomena, un- 
der a mathematical expression nowadays called Fluctuation Relation (FR). They 
proposed and tested the following formula 

P(A, 6 -As) 

where is the average of the power dissipated on evolution segments of dura- 
tion r, and PiA^ e ±Ag) is the probability that A^ takes a value in an interval 
of length S, centered at ±A. In analogy with the periodic orbit expansions [3,4], 
Eq.(l.l), was obtained from the "Lyapunov weights" in the long r limit. Remark- 
ably, Eq. (1.1), does not contain any adjustable parameter. 

In 1994, Evans and Searles obtained the first transient FRs for the energy 
dissipation rate divided by k^T, denoted by f2 [5, 6, 7, 8, 9, 10, 11], which we 
call transient Q-FR , and which concern the statistics of an evolving ensemble 
of systems, instead of the steady state statistics. The only requirement for the 
transient 's to hold is the reversibility of the microscopic dynamics. Because 
they describe the behaviour of O, these relations can be experimentally verified 
[12]. 

In 1995, Gallavotti and Cohen provided a mathematical justification of the 
Lyapunov weights of Ref.[2], introducing the Chaotic Hypothesis [13, 14, 15], 
which states that: "A reversible many-particle system in a stationary state can be 
regarded as a transitive Anosov system, ' for the purpose of computing its macro- 
scopic properties". The result was a genuine steady state FR, which we call A-FR , 
as it concerns the fluctuations of the phase space contraction rate A. This quantity 
equals the energy dissipation rate Q in a subclass of Gaussian isoenergetic parti- 
cle systems, which includes the model of [2], while in many other cases A differs 
from n by a total derivative. As far as we know, the works [13, 14, 15] provide 
the only answer that has been given so far to the question of which dynamical 
systems can be proven to verify the steady state A-FR. A strong assumption as 
the Chaotic Hypothesis raises the question of which systems of physical interest 
are "Anosov-like", since almost none of them is actually Anosov. The answer of 



Anosov systems are smooth and uniformly hyperbolic. 



Introduction 



3 



Refs.[13, 14] is that the Anosov property, in analogy with the Ergodic property, 

holds "in practice" for sufficiently chaotic, reversible systems. 

A complementary approach, which addresses a different question, has been 
recently proposed in [16], developing ideas first introduced by Evans and Searles, 
see e.g. [5]. The question concerns the physical mechanisms at work in systems 
which do obey the steady state Q-FR. This approach leads to a general frame- 
work within which various transient and steady state relations can be obtained 
for a variety of observables, and especially leads to the identification of physical 
mechanisms and time scales underlying the validity of the steady state FRs. The 
result is that, while the transient relations only rest on time reversibility, the steady 
state relations do require further properties to hold.^ In [16], a property called "Q- 
autocorrelation decay" is identified as the reason for the steady state Q-¥R , as 
well as of a number of other asymptotic relations, to be verified, when they do. 
This approach leads to the identification of the time scales concerning the FRs 
with the decay times of a quantity which we introduce in Section 3.3. Therefore, 
the need arises to test on concrete systems the validity of this identification and to 
quantify the corresponding time scales.^ These issues are closely related with the 
role played in the dynamics by the singularities of the total derivative which distin- 
guishes A from Q. This has been discussed in various papers, like [23, 17, 24, 19]. 
The present paper is devoted to the investigation of such questions. 

We study the heat transport and its transient and steady state fluctuations, in a 
mechanical dissipative system that is maintained out of equilibrium by an imposed 
temperature gradient. This local thermostating mechanism leaves unaltered the 
Newtonian dynamics in the bulk of the system, dissipating only at its boundaries. 
However, due to the interactions, dissipation may depend on the non thermostated 
degrees of freedom (d.o.f.) as well [22]. To impose a temperature gradient we use 
deterministic Nose-Hoover thermostats, that add a time-dependent friction term 
to the Hamiltonian of the system in such a way that, if the system is ergodic, the 
energy distribution of the thermostated d.o.f. is Boltzmann distributed with well 
defined temperature. As a consequence of the use of Nose-Hoover thermostats, 
the phase space contraction rate is unbounded, possibly leading to statistical prop- 
erties of its large fluctuations which are not described by the standard FR (1.1). 
This scenario has been considered in various works, and depends on the form of 
the tails of the distribution of the observables of interest, see e.g. [17, 19, 23, 24]. 

^For instance, the systems under consideration must converge towards a steady state. 

similar study of the time scales was performed in Ref.[18], in one investigation on the 
applicability of the Chaotic Hypothesis. 
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As far as particle systems are concerned, the corresponding modified FR appears 
to be consistent with the Nose-Hoover thermostated Lorentz gas [25]. In our case, 
consistently with the cited literature, the unboundedness of A does not imply any 
modification of the standard FR. The distribution of the fluctuations of A, indeed, 
appear to have Gaussian tails. Therefore, the search for deterministic particle sys- 
tems which verify a modified version of FR, apart from the case of [25], must 
continue. 

The paper is organized as follows: after defining the many-body dissipative 
model that we will use, in Section 2 we discuss the equilibrium and nonequilib- 
rium states that arise from the coupling with the local thermostats. In Section 3 
we study the transient Q-FR and the steady state fluctuation relations for the phase 
space contraction rate, and for the dissipation function. Our conclusions are pre- 
sented in Section 4. 



The model we consider consists of a gas of non-interacting point-like particles of 
mass m that move freely inside a channel made of L identical two-dimensional 
cells. Each cell consists of a rectangular region of height Ay and width Ax, 
containing two fixed freely rotating disks of radius R and moment of inertia 6. 

The position of the disks in each cell and their radius are chosen so that the 
geometry of the channel corresponds to a periodic Lorentz gas in a triangular 
lattice at critical horizon, defined as the smallest disk radius to disks separation 
ratio such that the length of the particle's trajectory between two collisions with 
the disks is bounded. The lattice of scatterers is shown in Fig. 1. Moreover, 
the particles and disks are allowed to exchange energy at collisions, according to 
"perfectly rough" collisions, that are reversible, conserve total energy and angular 
momentum. These collisions are defined by the following formulae, which relate 
the normal and tangential components of the particle's velocity v with respect to 
the disk's surface, and the disk's angular velocity u before (unprimed quantities) 
and after (primed quantities) the collision [26] 



2 The Model 



- V, 



n 




(2.1) 



Rjj' 



2 



1 + ?7 
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Thus, collisions are deterministic, time-reversible and phase space volume pre- 
serving. The parameter 77, defined as 

is the only relevant dimensionless parameter characterizing the collision. It de- 
termines the energy transfer between disks and particles in a collision. For finite 
values of 77, the particles indirectly exchange energy with each other, through the 
interactions with the disks, even though they do not directly interact. 

This model, introduced in [26] as a mechanical model for coupled heat and 
matter transport, shows realistic equilibrium and out of equilibrium thermody- 
namical properties, in very different situations, e.g. when subjected to external 
gradients of temperature and chemical potential, in the presence of external elec- 
tric or magnetic fields [27]. 

The simple energy exchange mechanism allows the equipartition of the energy 
among all the d.o.f., allowing the system to reach a local thermodynamical equi- 
librium. In particular, in microcanonical simulations, the particle velocities are 
Maxwellian distributed at uniform "temperature" consistent with the equipartition 
theorem. These temperatures also characterize the distribution of angular veloc- 
ities of the rotating scatterers. The same is true in canonical and grancanonical 
simulations, where the system is subjected to temperature and chemical poten- 
tial gradients. One consequence of the establishment of local thermal equilibrium 
is that the gas of particles is described locally by the equations of state of ideal 
gases. Recently, the properties of the stationary state for a modified version of 
this model, coupled to stochastic energy and particles reservoirs, have been ana- 
lytically obtained in the weak coupling limit [28], and finite coupling corrections 
have been obtained in [29]. 

2.1 Local Nose-Hoover thermostats 

To establish a nonequilibrium state, two heat baths are placed at the left and right 
boundaries of the model, with local temperatures set to Tq and respectively. In 
previous works [26, 27, 28, 29], the bulk of the system was coupled to stochastic 
thermochemical baths, so to exchange energy and particles at the boundaries. Dif- 
ferently, we consider a closed system (particle number is conserved), and model 
the heat baths as local Nose-Hoover thermostats [30]. Roughly speaking a Nose- 
Hoover thermostat adds a kind of friction term to a number of degrees of freedom, 
whose energy fluctuations turn to be canonical at some temperature T, as a result. 
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Figure 1: Geometry of a system of length L = 7 cells. The disks of radius R 
form a triangular lattice at critical horizon, meaning that the height and width of 
the cells equal, respectively, Aa; = 2R and Ay = 2W, where W = AR/^/S 
is the separation between the centers of the disks. Periodic boundary conditions 
along the vertical direction are considered. The dashed lines indicate the single 
hexagonal cells. The solid lines are hard walls. At the left top we schematically 
show a single collision: when a particle collides with a disk at an angle a with 
respect to the normal at the collision point, the outgoing angle of its trajectory a' 
is determined by Eq. 2.1. In between collisions, the angular velocity of the disk to 
evolves freely for the disks in the bulk (empty disks) and according to Eq. 2.5 for 
the thermostated disks (filled disks). 



More precisely, consider d.o.f. with generalized positions q, momenta vr 
and mass m, subjected to an external potential 0(g) described by the Nose Hamil- 
tonian [31], 

ff, Jg, vr, X, S, s) = 5^ + 0(g-) + 

S (2.3) 

+ =^ + {g + l)hT\nX . 

Here X and S are the generalized position and momentum of the thermostat and Q 
its effective mass, s is a generalized time variable, T the temperature, g a number 
related to n,u and k^, the Boltzmann constant. If one identifies s and vr with the 
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physical time t and momenta p, letting ds = Xdt and tt^ = Xpj, and additionally 
g = then, it is easy to show that if the generalized variables tt, X and S are 
distributed microcanonically, then the physical variables q and p are distributed 
according to the canonical distribution with temperature T [32]. 

Furthermore, if one is solely interested in real time averages, the equations of 
motion that derive from (2.3) can be rewritten in terms of the physical variables q, 
p and t, eliminating the variables X, S, tt and s to obtain [30], 



(2.4) 

th 2 ^ 




where r^j is the relaxation time of the thermostat related to its effective mass Q 
as — Q/n^i^hgT and ^ — E/Q, (in what follows we will refer to the variable 
^ as the thermostat). Inspecting Eqs. (2.4), known as the Nose-Hoover equations 
of motion, the role of the thermostat ^ becomes evident: it acts as a friction term 
that appropriately changes sign, depending on whether the instantaneous kinetic 
energy of the thermostated d.o.f. is larger or smaller than the desired mean 
kinetic energy 

To simulate a nonequilibrium state, we couple Nose-Hoover thermostats to the 
angular velocity of each of the four disks at the boundaries of the channel. Since 
the disks are pinned to their positions, the thermostats can interact with the rest 
of the d.o.f. only via the collisions with the point-like moving gas particles. We 
denote by and the thermostats for the upper and lower disks at the cold side 
and as E,^ and the thermostats for the upper and lower disks at the hot side, thus 

= 1. Note that it is possible to couple the two disks at each boundary to one 
single thermostat, in which case, = 2. We have chosen the former coupling as 
we get more efficient numerics than in the latter case. Also, we assume that there 
is no external potential acting on the particles. 

In [33], it has been shown that the Nose-Hoover dynamics are not ergodic for 
= 0, thus the phase space probability distribution is not canonical. However, 
as we will see, the interaction of the thermostated disks with the gas of particles 
leads to ergodic motions in the accessible phase space. 

Labeling the cells of the channel from 1 to L, the equations of motion for our 
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model can be written as follows 



Pi 



Pi 
m 



l,...,n 



1, 
2, 



UJr 



± 



± 



± 







'th 




^th 





- 1 



(2.5) 



4^ 



where n is the number of particles. The abstract operators T^{t), represent the 
instantaneous forces exerted on particles and disks during a collision. They are 
subjected to the condition 



i=l 



E ^^r,;(t) + E ®''^^V■,^(t) + E ® (^cT,;c(^) + ^hT,;h(^)) = , (2.6) 



L-1 



j=2 
{+-} 



{+,-} 



which accounts for the conservation of energy during the collisions. 

The dynamics (2.5) is dissipative as the phase space volume contracts at a rate 
given by 

A = -divpf = + + + , (2.7) 

where V = (q, p, w, ^). The mean value of the phase space contraction rate A is 
zero only at equilibrium, i.e., when = T^. As we will see in Section 3.2, 
the cold thermostats contract the phase space volume at a rate faster than the 
expansion produced by the hot thermostats so that on average, A is positive. 

2.2 Ergodicity and Equilibrium State 

Because our system is not subjected to any external field, the thermostats dissipate 
energy only when they collide with the particles. These collisions, in turn, make 
apparently ergodic the evolution of the thermostated d.o.f., since the particles dy- 
namics is randomized by the motion in the bulk of the system. To illustrate this, 
we have considered the dynamics inside one single hexagonal cell (see Fig 1), 
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containing a thermostated disk and n = 10 particles, with reflecting boundary 
conditions. In Fig. 2, the evolution of the thermostat ^ and the angular veloc- 
ity of the disk uj is shown for a temperature T = 1000. In the main panel, the 
solid curves correspond to the evolution of the isolated thermostated disk, i.e., 
without the interaction with the gas of particles, for three pairs of symmetric ini- 
tial conditions. In this case, the evolution of the vector {u,^) is oscillatory: u 
oscillates around its equilibrium value ^/T, preserving its initial sign. The ther- 
mostat ^ oscillates symmetrically around zero. Differently, when particles are 
allowed to interact with the disk, the trajectory of any initial condition explores 
all the available phase space, as shown by the cloud of points in the main panel 
of Fig. 2. Furthermore, the histograms in the left and lower panels correspond to 
the marginal distributions of the points of the trajectory, which coincide with the 
corresponding equilibrium distributions of ^ and uj (dashed curves). 

Summarizing, while in the absence of interactions the dynamics of the disk 
is trivial (periodic), the interactions with the particles lead to trajectories which 
appear to explore all the available {u, space. This apparent ergodicity, combined 
with the observation that the distribution of the particles velocities approximates 
a Maxwellian at the temperature T, indicates that locally (in the cell) a form of 
thermal equilibrium is established. 

In the single cell we also have computed the energy probability distribution 
function P{E) of the disk and of the energy of the gas of particles, for a given 
temperature T = 100. The results are shown in Fig. 3. The agreement with 
the theoretical equilibrium Boltzmann distributions is excellent. Independently of 
the initial state, the gas of particles and the thermostated disk reach quite rapidly 
(OiW^) collisions times) the equilibrium state at a temperature equal to that of 
the thermostat. For the model depicted in Fig. 1 a uniform equilibrium state is 
estabUshed if Tc = T^ = T. 

A similarly fast relaxation to equilibrium, due to the interaction of Eq. (2.1), 
was also observed using stochastic baths [26, 27]. The local equilibration between 
the thermostated and not thermostated d.o.f. is due to the fact that, at the collisions, 
the incoming and outgoing temperatures dynamically relax to the same value. 
In the inset of Fig. 3 the incoming and outgoing energy probability distribution 
function conditioned to the occurrence of a collision is shown. Both PDF coincide. 
Moreover, they coincide with the theoretical equilibrium conditional probability 
(shown as a solid line). This is a remarkable feature of the present model which is 
not always obtained in other systems [34]. 

Note that the same qualitative results are found independently of the particle's 
density n. However, the time required to reach equilibrium with the gas of parti- 
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Figure 2: Main panel: Evolution of the angular velocity of the disk uj{t) and its 
thermostat ^(t) for the one single cell at a temperature T = 1000, in the pres- 
ence (dots) and absence (solid curves) of the gas of particles. In the presence 
of the particles the trajectory is obtained from one single initial condition. The 
solid curves correspond to three pairs of symmetric initial conditions. Secondary 
panels: histograms of the marginal distributions P(Q (left) and Pico) (lower) of 
the trajectory with interactions. The dashed curves represent the corresponding 
equilibrium distributions. 



cles is determined by the rate 7 oc n^/T at which the gas of particles collides with 
the thermostated disk, compared to the relaxation time % of the thermostat. 

Also note that the thermostated d.o.f. to changes discontinuously at collision, 
while the thermostat variable ^ does not. We have not observed any noticeable 
consequence of the uncoupled evolution of ^ and to at the collisions. We expect 
that this will affect the efficiency of the thermostat, at most. A similar observation 
was made in [33], for a thermostat with same kind of discontinuities. 
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Figure 3: Equilibrium energy probability distribution P(E) for the disk (crosses) 
and particles (circles). The distributions were obtained for 10 particles in a hexag- 
onal cell containing one single disk thermostated at a temperature T = 100. The 
dashed lines correspond to the theoretical equilibrium distributions. Inset: Condi- 
tional probability distribution of the energy of the particles P{E^, before (circles) 
and after (plus) the collisions. The dashed line corresponds to the theoretical dis- 
tribution function. 



2.3 Nonequilibrium Steady State 

We now turn our attention to nonequilibrium states. In [26, 27] the transport prop- 
erties of the system of Sec. 2 were studied. The system was coupled to stochastic 
thermochemical baths to create temperature and chemical potential gradients. In 
that configuration, it was found that the hypothesis of local thermal equilibrium 
holds and that the deviations from it are consistent with the transport of heat and 
matter that appears as a response to the external gradients. LTE results from the 
mixing properties of the effective interaction among particles. In the absence of 
interaction, e.g., when ?7 = 0, LTE is not satisfied (see e.^.,[35]). 

In our present situation, the stochastic baths are substituted by deterministic 
Nose-Hoover thermostats. Moreover, the Lorentz channel is closed, i.e., there is 
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Figure 4: Probability distribution function of the angular velocity P(uj) of the 
"cold" (circles) and hot (squares) thermostated disks, from a simulation with Tq = 
900 and T^j = 1100, n = 20 particles per cell and L = 7 cells. The dashed lines 
correspond to the theoretical equilibrium distribution functions. In the inset, the 
same distributions are shown, for a much larger temperature gradient of Tq = 100 
(circles) and Tjj = 1900. 



no exchange of particles between the system and the reservoirs, implying that the 
particle current is zero. 

Since LTE is a consequence of the mixing due to the interactions, we have 
observed that it holds, indeed, in the bulk."^ However, away from equilibrium, 
this is not necessarily the case in the thermostated regions. In order to understand 
how effectively the interactions establish LTE in the neighbourhood of the ther- 
mostated disks, we have performed extensive numerical simulations, at different 
temperature gradients. In Fig. 4 we show the probability distribution function of 
the angular momentum of the thermostated disks P(uj), for a temperature gradi- 
ent of = 900 and = 1100. The circles and the squares correspond to the 
distribution of the disks at the cold and hot sides, respectively. The dashed lines 

'^By bulk we mean all cells which do not contain a thermostated disk (open disks in Fig. 1. 
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correspond to the theoretical equilibrium distributions at the nominal tempera- 
tures of the thermostats and Tj^. The agreement of P(uj) with their respective 
theoretical distributions is excellent. 

However, for very large temperature gradients, LTE is not established at the 
cold side, while the LTE at the hot side does not necessarily correspond to the 
chosen T^. For instance, the inset of Fig. 4 shows the cold and hot P(uj) for Tq — 
100 and = 1900. The exponential tails of P(uj) at the cold side is a signature 
of the bad equilibration. While at the hot side P{uj) remains Gaussian, but its 
variance corresponds to a temperature which is not T^. The bad equilibration 
around the cold side is due to the inefficiency of the thermostat to thermalize the 
much more energetic particles coming from the hot region. Thus, the bulk of 
the system practically equilibrates with the hot thermostat alone. This problem 
can be partially solved by varying the relaxation time of the cold thermostat Tji,. 
However, the state at the cold side is very sensitive to variations of which 
must be determined case by case. In the following we restrict our computations to 
gradients AT/T < 1, for which LTE, at the correct temperatures is established. 

In all cases we have found that the steady state is characterized by linear par- 
ticle n(x) and energy T{x) density profiles. In Fig. 5, we show an example of 
n(x) and T(x) for a chain of 48 cells and 10 particles per cell, with Tq = 800 
and Tfj — 1200. Note that the temperature gradient generates a particle density 
gradient, in the bulk. However, since the billiard model that we consider is closed, 
the mean particle current is zero. Note also that T(x) does not jump at the bound- 
aries, which means that a good equilibration between thermostats and gas particles 
is achieved. By measuring the heat current as a function of the size of the system, 
for a fixed temperature difference, we verified the validity of Fourier's law to high 
accuracy. 

The same results had been previously obtained in [26, 27], for the system cou- 
pled to stochastic thermochemical baths. It is indeed expected that the transport 
properties in the bulk of the system do not depend on the particular model of heat 
bath. 

3 Fluctuation Relation for Heat Conduction 

In this section we study the fluctuations of the phase space contraction rate and of 
the entropy production, or energy dissipation in the transient response to external 
drivings and in the stationary state of the model system introduced in Sec. 2. 
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Figure 5: Stationary profiles of the density n{x) and the energy density T(x) of 
the particles from a simulation with = 800 and = 1200, n = 10 particles 
per cell and L = 48 cells. The dashed lines correspond to a linear profile between 
T(j and Tj^, (lower panel), and to a best linear fit of n{x), (upper panel). Note, 
however, that for our system neither T(x) nor n(x) need to be linear. 



3.1 Transient Fluctuation Relation 

Consider a system that is initially at equilibrium. At time t = 0, an external force, 
either mechanical or thermodynamical, is applied. For times t > the system 
evolves from its equilibrium state towards, possibly, a nonequilibrium steady state. 
In the process, the system dissipates the work that the external forcing performs 
on it. Furthermore, consider an ensemble of trajectories in the phase space F, that 
originate at the equilibrium state (at t = 0) and evolve for a finite time r. For 
S > 0, let Aj = (A~ S,A + S) and Aj = (-A - 5,A + S). The Transient 
Fluctuation Relation [36], states that 

1 Prob {E{n,^, € Aj)) 

m -. — =^ f ~ , (3.1) 

Prob (E(fio,. e AJ)) 
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where E(0,q ^ G (a, b)) is the set of points in F, such that Qq ,^ G (a, b) ^ and 
r^Q^ = (l/r) Jq^ ri(F(t))(it is the time averaged dissipation function, which is 
defined in terms of the (initial) phase space probability density / and of the phase 
space contraction rate A as [36]: 

no,,(F) = ^ (^In + A(m)d?j . (3.2) 

Remark 3.1 The density f can be arbitrarily chosen. However, if it represents 
the equilibrium ensemble associated with the chosen dynamics, VL corresponds to 
the physical dissipation (the entropy production, close to equilibrium) [16]. 

In this section we study the fluctuations of the time averaged dissipation func- 
tion as the system moves away from a given equilibrium state. In Ref.[37], a tran- 
sient fluctuation formula for heat conduction was obtained, considering a generic 
system coupled to two (deterministic) thermostats, which is the case of the present 
work. In Ref.[16] the transient fluctuation relation has then been obtained in great 
generality for under the sole assumption of time reversibility. Hence we recall 
the definitions of [16], which will be used for the steady state fluctuation relations 
as well. 

Consider a dynamical system with phase-space M, described by an evolution 
operator 5* : — > 7W. Let S* be reversible, i.e. such that for an involution 
i : Ai ^ M. representing the time inversion operator, iS^T = S^^iT holds for 
all F G and all t G M. Introduce a probability measure fi with density / on A^, 
i.e. d/x(F) = f(T)dT, which is even for i, i.e. f(iT) = /(F). Let : M^Rhe 
any odd observable, i.e., 0(iF) = —0(F), of a many-body dynamical system that 
satisfies the assumptions above and denote by 

^0,r = l r<P(S'T)dt, 







the time average of during a time r with initial condition F. 

In [16] it was proved that (f) satisfies the following symmetry relation: 

where f2Q^(F) = rno,r(r) is the time integral of over the interval [0,r]. In 
the r.h.s. of (3.3) the average (■)^^ is the ensemble average over the set of 



In what follows, we will omit the specification of the set E for the sake of simplicity. 
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trajectories that satisfy the constraint that 0o,T ^ ^5^- The Conditional reversibility 
theorem derived in Ref.[38] is related to the steady state version of (3.3), given in 
[16]. 

When the system is subjected to an external mechanical force F, the density 
function leading to the physical dissipation is the equilibrium {F — 0) probability 
density. Other choices are possible and, for instance, the uniform density in a 
compact phase space, /(F) = 1/|A^|, yields Vt — K, the phase space contraction 
rate. 

Remark 3.2 The relation (3.3) is exact. It is purely a consequence of the time 
reversibility of the dynamics. Remarkably, this relation is valid for any observable 
that is odd with respect to time reversal. 

In particular, we can take I^q^ as the observable. Doing so, Eq. (3.3) becomes 



where the error term e < S, in general will depend on 5, A and r. In [16] the 
transient fluctuation relation (3.4) was called il-FR . 

3.2 Dissipation Function for Heat Flow 

In this section we investigate the transient response of the system of Sec. 2 to a 
temperature gradient. The dynamics of the model is given by the equations of 
motion (2.5) and the energy density of the system H, can be written as 

1 " f) ^ 

= i E + ^ E + (-7)'] + 

i=l j=l 



{+,-} 



(3.5) 



where uf = and uj^ = uj^. Here, = + ^nd = + ^h- We 
consider the following transient process: 

At times t < the system and the thermostats are in thermal equilibrium at 
some temperature Tq. For simplicity we choose the equilibrium temperature as 
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Tq = (Tq + Th)/2 (in the linear regime Tq is the mean temperature of the steady 
state of the system). Therefore, the initial density distribution corresponds to the 
canonical equilibrium distribution 



with /?o = I/^bTq, and H(r) the energy density (3.5), with = = Tq. At 
time t = the temperatures of the thermostats are set to for the cold thermostat 
and Tjj for the hot one. At times i > the system is not in equilibrium and a heat 
flux develops. 

We measure the transient response of the system by considering time averages 
during a finite time r of the dissipation function Q. 

Substituting (3.6) into (3.2), the time average of fig.r takes the simple form 

^o,r(r) = ^ {PoH(S'^) + MS'r)) dt = ( HiS'D-HiD ^^ ^ ^^^^^^^ 

(3.7) 

The quantities and A, thus differ by a total time derivative, as in the cases 
considered in e.g. [17, 18, 24, 25]. Taking the time derivative of (3.5) and using 
the constraint (2.6) we obtain 

Hit) = - (TcCc + ThCh) ■ (3.8) 

Finally, using the expression for the phase space contraction rate (2.7), the 
time averaged dissipation function for the heat flow of our model becomes 

no,,(r) = (Cc; 0,. - 0,.) ■ (3.9) 

As a test of the physical character of the dissipation function (3.9), we have 
computed its long time average in the steady state, Qo,t' a function of the rela- 
tive temperature gradient /ST /T. This is shown in Fig. 6 for a chain of 5 cells and 
20 particles per cell. For each simulation the temperatures of the thermostats Tq 
and Tjj were varied symmetrically around T = 1000. At equilibrium fip r = 0- 
Out of equilibrium, grows with the external gradient as (AT/T)^. In the inset 
of Fig. 6 we show the average behaviour of the individual thermostats and ^y. 
At equilibrium, both thermostats average to zero. Away from equilibrium, their 
mean values grow in absolute value, with the distance from equilibrium. While 
the dissipation at the cold thermostat contracts the phase space volumes, the hot 
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Figure 6: Long time average in the steady state of the dissipation function 
= ^(^c ~ ^h) ^ function of the relative temperature gradient AT/T. The 
dashed line corresponds to (AT/T)^. Inset: local phase space contraction rate, 
(circles), and (stars). 



thermostat expands them. However, dissipates more than creates, leading, 
on average, to a positive global dissipation in the steady state. 

To test the Transient Fluctuation Relation (3.4), we have numerically com- 
puted the probability distribution function of the dissipation function (3.9), for 
different values of the transient time r, and for a channel of L = 10 cells with 
n = 4 particles per cell. This test may seem unnecessary, because our dynamics 
is time reversal invariant, which is the only condition required for (3.4) to hold. In 
reality, there are two subtle points, at least, which make non-trivial the test: a) our 
simulations are not exactly reversible, because of the numerical scheme approxi- 
mating the solution of the equations of motion, and because of round-off errors; 
b) the initial continuous ensemble /, required by the theory, is replaced by a finite 
number of initial states, supposed to be picked up at random with density /. This, 
combined with the fact that the dissipation takes place discontinuously in time, 
because it relies on the interactions between thermostated disks and point parti- 
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cles, makes delicate the choice of the numerical procedure, which could introduce 

spurious effects, leading to fluctuations that do not satisfy the Transient Fluctua- 
tion Relation. Such spurious effects could spoil also the tests of the Steady State 
Fluctuation Relation which, differently from the Transient Fluctuation Relation, 
requires more than time reversibility to hold, and certainly needs to be tested case 
by case. 

Therefore, the test of the Transient Fluctuation Relation is rather important, in 
our system, to calibrate the numerical protocols. Our investigation has shown that 
the effect of the numerical integration on time reversibility is negligible, as far 
as (3.4) is concerned, but also that the construction of the initial ensemble needs 
special care. What proved adequate to our purpose is the construction of / from an 
equilibrium simulation, which generates a collection of decorrelated microstates, 
by saving a microstate every sufficiently long time. 

We have proceed as follows: at t = the state of the system is set to a 
microstate previously generated by the equilibrium dynamics at a temperature 
Tq — (Tq + Th)/2. From this initial state, the temperatures of the thermostats 
are set to Tq and respectively. Then, the system evolved for a certain time 
r, and CIq ^ was computed. Considering a large ensemble of such processes, the 
empirical probability distribution function P(r2Q^) was constructed. 

Furthermore, we define 

Ca(r; A) = -\n 4 , (3-10) 

where A is the centre of the corresponding bin of the empirical distribution. In 
terms of (3.10), the transient Q-FR is written as 

C^(t;A) = A; for all r. (3.11) 

We have verified the Transient Huctuation Relation for several values of r e 
[0.01, 20], for a temperature gradient given by = 700 and = 1300. These 
results are shown in Fig. 7 for the slope of the relation (3.11). In the inset of 
Fig. 7 we show P(ilQ .^) for four values of r. We have found that for small r, the 
obtained distribution is highly asymmetric. It becomes Gaussian only for large r. 

The verification of (3. 1 1) indicates that our numerical procedure is appropriate 
to test the 0-FRs, hence that it may be used to test the Steady State Fluctuation 
Relation. 




Figure 7: The slope of Cq(A, t) as a function of r, for a channel of 10 cells, 
Tc = 700 and = 1300. Inset: Probability distribution function P(Uo r) for 
(from left to right) r = 1, 5, 10 and 15. The dashed line corresponds to a Gaussian 
function. 



3.3 Steady State Fluctuation Relation 

The steady state fluctuation relation obtained by Evans, Cohen and Morris s in [2] 
was motivated by the fact that the phase space contraction rate A is proportional to 
the energy dissipation rate divided by the temperature, for a subclass of Gaussian 
isoenergetic particle systems, which includes the model of [2]. This way the fluc- 
tuation relation for A, which we call A-FR as in [16], could be given a physical 
meaning in [2]. The A-FR has been later proved by Gallavotti and Cohen to hold 
for all dissipative, time reversible, smooth and uniformly hyperbolic dynamical 
systems, independently of their physical meaning [13, 14]. For such systems, the 
A-FR states that: 

lim log = A , (3.12) 
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where = ^o,t/ (^)oo phase space contraction rate averaged during a time 
interval r and normalized to its stationary state value (A)^^, and P(A) is the cor- 
responding steady state probability distribution function. The A-FR is supposed 
to be valid only if (A)^ 7^ and A belongs to the interval of definition of (A)^^. 

In [16], it has been shown that the steady state fluctuation relations may only 
rely on the time reversibility of the microscopic dynamics and on the time be- 
haviour of the quantity B^, defined in Eq. (3. 14) below, independently of the fluc- 
tuating odd observable under consideration. This way, O appears to be a sort of 
"thermodynamic potential", which determines the behaviour of all other odd ob- 
servables, but works even far from equilibrium.^ This can be understood starting 
from the symmetry relation (3.3), concerning the fluctuations of any observable 
that is odd with respect to time reversal. In (3.3), the symmetry of the fluctuations 
is expressed in terms of averages of the dissipation function (3.2), with respect to 
the initial measure /i of density /. Taking the ensemble average with respect to the 
measure evolved up to time t, n^, of the averages of il along trajectory segments 
of duration r, which at time t — started from the initial state one obtains: 

1 /i, (Hn , G At) 1 

- In = A + 6(5, t, A,T)-- InBait, r) , (3.13) 

where e is bounded hy < e < \5\, and B^it, r) is a kind of correlation function 
given by^ 

B^(t, T) = (e-^o.e-".+..*+.)^^^^^^^^+ . (3.14) 

If fj,^ tends to a steady state fj,^ when t — > 00, then (3.13) should change from a 
statement on the ensemble at time t, to a statement on the statistics generated by a 
single typical trajectory of the stationary state This requires some assumption, 
because t tends to infinity before r does and, in principle, the growth of t could 
make the conditional average in (3.14) diverge. However, if a finite time scale t„ 

^In the mathematical theory of smooth, uniformly hyperbolic dynamical systems, it is A that 
plays a role analogue to that of a "thermodynamic potential". 

''Although the aspect of ^^(t, r) resembles that of a standard correlation function, it cannot 
be always interpreted like that. For instance, in cases which verify the modified FR of axiom C 
systems, which enjoy a strong form of correlations decay [40, 39, 41, 42], the term In BqH, t) 
does not decay, if 17 is chosen to equal A. In other cases in which no FR is verified, such as those 
with trivial steady states considered in [21, 20], BQit,T) grows exponentially with t, as simple 
calculations show. Because the t ^ 00 Umit must be taken before the r ^ 00 Umit, this illustrates 
that the transient FR cannot not turn into a steady state FR, in these cases, although the transient 
FR holds at all times. A more detailed discussion of these facts can be found in [1, 16]. 
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characterizes the decay of In Bq,^ then (3.13) shows that is the time scale 
which must be exceeded for the t ^ oo limit to be practically reached. When 
this is the case, the transient relations yield the steady state ones, in the long time 
limits [16]. 

In some cases, the steady state A-FR and Q-FR coincide, as mentioned above, 
while in other cases they are clearly different. For instance, in some cases in which 

A represents the heat flow, Q the dissipated energy, and the particles interaction 
potentials are singular, the tails of the A-FR and Q-FR may be quite different 
[17, 16, 23, 24] ^. Which systems actually enjoy the decay of In B^it^r) 
required for (3.12) to hold, or for other steady state relations to be derived from 
(3.3) is not known. Hence a test is necessary in the systems of physical interest 
Therefore, we perform the test on our model. 

The mesoscopic time scale r„j, expresses a physical property of the system, 
typically the mean free time and thus, does not depend on t or r. If t, r 3> r^, the 
boundary terms f2j_^ ^ and , are typically small compared to ^t,t+T- 

For this not to be the case, some singularity of Q must occur within (t — t^, t) or 
(t + T,t + T + r„). However, even in that case, similar events may equally likely 
occur in the intervals (0, t) and (t + r, 2t + r), hence ^ and are 
expected to contribute only a fraction of order 0(t^/t) to the arguments of the 
exponentials in (3.14). Therefore, in the cases of thermodynamic interest, one can 
write 



with an accuracy which improves when t and r are growing multiples of the time 
T„^. The last line is just a statement of the conservation of probability and of 
the fact that Bq typically decays within microscopic times [1]. One then realizes 
that B^it, t) tends to approximately 1 as 1/r, with a characteristic scale of order 

^Note that Bq is computed with respect to the initial measure fi, even for t > 
^This behaviour was first obtained for a stochastic model, namely the Brownian particle 
dragged in a liquid by a moving harmonic potential [23] . Later, the corresponding extended A-FR 
has been derived for observables for which the probability distribution function has exponential 
tails [24]. 

'"it is known, however, that r^^ In Bq has to decay in systems not too far from equilibrium, i.e. 
in the linear regime of Irreversible Thermodynamics, for the transport coefficients to exist [1, 16]. 





'o,t+Tm \ /e.~^t+T+Tm,2t+T 



(3.15) 



0(1), 
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0{t^), which leads to the validity of the steady state Q-FR . Then, for sufficiently 
small 5, and large t and r the relation 

iln^^42^^^±^^^«A. (3.16) 

holds, and remains valid in the t ^ oo limit. If P is the corresponding steady 
state probability distribution function, one can rewrite (3.16) as 

1 P (Hn , e At) 

Ca (r; A) = -ln )_ ^ « A ; for r > , (3.17) 

which we call steady state Q-FR. 

We have measured the steady state fluctuations of A and of O, for the model 
introduced in Section 2. Let us first describe those of A, and then compare them 
with those of ^2, in the r ^ oo limit. 

We consider a channel of length L = 10 cells and n = 4 particles per cell. As 
before, the disks in the leftmost and rightmost cells are coupled to Nose-Hoover 
thermostats with = 700 and = 1300 respectively, and T^^^ — 1. Starting 
from a random initial state we leave the system relax for a time corresponding 
to 5 X 10^ collisions of the particles with the disks. This relaxation time ensures 
that the system has reached its stationary state at which, the mean heat current 
is uniform. Using Eq. (2.7), in the stationary state we measure the values of the 
phase space contraction rate A each 10~^ time units. We then compute the time 
average Aq^^, from disconnected intervals of the time series of A. Finally, we 
compute the empirical probability distribution function of Aq^, normalized to its 
steady state value (A)^. 

We have found that, for our system, the unboundness of A is not associated 
with exponential tails in its probability distribution. This is shown in Fig. 8-a, 
where P(A^) is plotted for different values of r: within numerical accuracy, P(A^) 
has Gaussian tails for all values of r. Note, however, that for the values of r 
considered, P(A^) is not a Gaussian distribution. 

To discuss the behaviour of the A-FR (3.12) for finite averaging times r, we 
define, as before, 

C^(r;A) = -i-ln^|^^^, (3.18) 

so that the relation (3.12)), for finite but sufficiently large r, would imply that 
Cj^(r; A)/ A 1. In panel h of Fig. 8, the dependence of Ca(t; ^4) on A is shown 
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Figure 8: Panel a: Probability distribution P(A^) for r = 1, 10, 20 and 40. The 
largest the r the more narrow the distribution gets. The data was obtained for a 
channel of L = 10 cells with a density of n = 4 particles per cell. The tempera- 
tures of the thermostats were set to = 700 and = 1300. Panel b: Cj^ir; A) 
as a function of A for different values of r: 0.1 (stars), 0.5 (squares), 1 (crosses), 
5 (circles) and 10 (pluses). The dashed line corresponds to the line with slope 1. 
As observed, C^{r; A) is approximately linear with A, for all r, with a slope ^ 1, 
that depends on r (see Fig. 9 below). Note that P(A^) are not symmetric about 
their maximum. Therefore, for the values of r considered, P(A^) is not Gaussian, 
but their tails decay as fast as Gaussian. 



for different values of r. Independently of the value of r, C^(t; A) is approxi- 
mately linear in A. The only dependence on r is the slope of C^(t; A) and the 
value of A at which the statistical errors become large. This is different from 
recent findings for other Nose-Hoover thermostated systems. In [25] the Nose- 
Hoover thermostated Lorentz gas was studied, and it was found that the fluctua- 
tions of A obeying the A-FR are only those that are smaller in magnitude than the 
steady state mean (A)^^^. Larger heat fluctuations follow the extended version of 
the A-FR of [23, 24], which saturates at some value. In these works, the deviation 
from the standard A-FR was supposed to be due to the singularities of A, which 
may produce exponential tails in the probability distribution function of A, pos- 
sibly as a consequence of the exponential tails of the distribution function of the 
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Figure 9: Comparison between the slopes of CJr; A) (squares), and Co (t; A) 
(circles), as a function of the averaging time r. The parameters are as indicated 
in Fig. 8. In the inset, Cj<^{t; A) /A is plotted in inverse units of r. The dashed line 
corresponds to the function e~^^/^. 



values of H. 

Our results confirm that the identification of deterministic particle systems 
which verify a modified version of the FR is far from obvious [17, 16]. Indeed, 
the most common situation which is found in the specialized literature, as well as 
here, is that even systems with singular A verify the standard FR, since the tails of 
the distribution of the fluctuations of A decay as fast as Gaussian tails. Then, as 
explained e.g. in Ref.[24], the standard FR must hold. In our case, the phase space 
contraction rate A is unbounded because of the logarithmic potential present in the 
Nose Hamiltonian (2.3), but the tails of the distributions of its fluctuations appear 
to be Gaussian, and C^(r; A) turns out to be linear in A, even for fluctuations much 
larger than (A)^. 

We now discuss the behaviour of Cj^(t] A) with r. In Fig. 9, the slope of the 
C/^{t; A) (squares), is shown as a function of r. These values were obtained from a 
best fit to a linear function of Cj^ir; A). For small times r, shorter than the meso- 
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scopic time scale r^, the slope of Cj^(t;A) decreases approximately as ~ l/r. 
The slope reaches a minimum value, after which Cj^ir; A)/A increases, seem- 
ingly converging to its asymptotic value 1. This convergence becomes clearer in 
the inset of Fig. 9, where the same data has been plotted as a function of l/r. 
These results are consistent with an exponential convergence of Cj^(t; A)/ A to 1 
as T — > oo (the dashed line in the inset of Fig. 9, corresponds to an exponential 
scaling of exp(— 26 / r)). Recalling Eq.(3.8), and the fast convergence of the 17 -FR, 
one may think that the monotonic convergence of the A-FR goes as l/r, rather 
than exponentially. Indeed, in the cases in which A is bounded, a conservative 
estimate of the convergence rate yields (Aj^ — A^)/r, where Ajy^ and A^ are the 
maximum and minimum values of A, see e.g. [17]. However, convergence may 
be faster even when A has no bounds, as in our case. Therefore, the steady state 
fluctuations of the phase space contraction rate for our model system satisfy the 
A-FR (3.12). 

Let us now turn to the steady state Q-FR (3.17) and compare its behaviour at 
large averaging times r with that of the fluctuations of the A-FR. Fig. 9 shows the 

slope of Co (r; A) (circles), as a function of r. For times shorter than the meso- 
scopic time scale r„j, the steady state fluctuation relations of A and of Q converge 
with an error which decays as l/r. At longer times (r > r^) the slope of Cq 
is fixed to 1, indicating that the steady state Q-FR (3.16) holds. Moreover, Fig.9 
shows that the convergence of Cq is much faster than that of C^, as expected 
because of the total derivative by which Q and A differ. Indeed, we observed that 
the turning point of the convergence of the A-FR occurs when Bq has decayed, 
and the Q-FR has converged, i.e. when only the term in H is left to decay. The 
difference in convergence rates is expected to grow when the system approaches 
equilibrium, and the dissipative part of A becomes dominated by the total deriva- 
tive ij." 

These observations support the idea of [16] that the scale set by the decay of 
IuBq determines the behaviour of the fluctuations of all odd observables, in 
our model, although the global behaviour of different observables may be differ- 
ent. To test this we have computed the behaviour of B^it, r). We considered a 

chain of L = 10 cells, density n = 4 particles per cell and Tq — 700, = 1300. 
As in Section 3.2, at t = the system is set to an equilibrium microstate at tem- 
perature Tq = (Tq + T^)/2. With the temperature of the thermostats set to Tq and 

"similar observations were made in [18], where the issue of convergence rates of FRs was 
thoroughly investigated in models of shearing fluids. In the case of [18], the convergence of the 
A-FR appeared so much slower than that of the fi-FR, that it could not be directly observed, and 
was only indirectly inferred. 
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Figure 10: Logarithm of B^it, r) of Eq. 3.14 as a function of the averaging time 
r, for t = 1 (solid squares) and t = 10 (open circles), for A = 2 and S = 0.1. The 
parameters are as indicated in Fig. 8. 



Tjj respectively, we follow the evolution of the system for a time 2t + r and com- 
pute the time averages of il^ ^, ^t,t+T ^t+T,2t+T- Finally, using Eq. (3.14) and 
a large ensemble of such processes, we obtain B^it, r). Fig. 10 shows the depen- 
dence of —T^^ \nBQ(t, t) on t for a = 2, 5 = 0.1, and t = 1 (solid squares) and 
t = 10 (open circles). For both values of t, —t~^ lni3f^(t, r) decays, and reaches 
zero within numerical accuracy, for times 10. This time scale coincides 

approximately with the time scale at which, /A converges to its asymptotic 
value 1, and Cf^/A starts its monotonic convergence. 

4 Conclusions 

We have discussed the behaviour of the transient and steady state fluctuations of 
the phase space contraction rate and of the dissipation function, in a Nose-Hoover 
thermostated system. The dynamics in the bulk of the system is purely Hamilto- 
nian, hence preserves the phase space volumes. Two Nose-Hoover thermostats, 
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are coupled to d.o.f. that are fixed at the boundaries of the system. 

We have studied the dynamics of the local thermostats and have characterized 
the equilibrium and nonequilibrium states of the system. In the equilibrium state, 
we have shown that the effective interaction among the particles is sufficient to 
ergodize the, otherwise regular, dynamics of the thermostat. This leads to a local 
thermalization of all the d.o.f. of the system, for /^T/T < 1 (which is our case). 
Out of equilibrium, when the temperatures of the thermostats are unequal, the 
system develops nonequilibrium temperature and density profiles. Consequently, 
a stationary and uniform heat current appears, transporting heat from the hot to 
the cold reservoir. Energy is then dissipated at a constant rate, measured by the 
dissipation function. 

We have found that the transient fluctuations of Vt do satisfy the transient Vt- 
FR, indicating that our numerical solution of the equations of motion and genera- 
tion of the equilibrium ensemble from the equilibrium dynamics are appropriate. 

For the system studied in this paper, in which A does not coincide with the 
energy dissipation rate Jl, the fluctuations of both quantities obey the same sym- 
metry given by the standard steady state FR, in the t ^ oo limit. In particular, the 
standard A-FR is verified despite the fact that A is unbounded, because the tails 
of the distribution of its fluctuations decay as fast as Gaussian tails. The A-FR 
needs In to have decayed, before starting a monotonic convergence, while 
the Q-FR appears to have converged, within numerical accuracy, after those times. 
This different behaviour depends on the total derivative H, which distinguishes A 
from Vl. The characteristic time scale of the decay of In is the mesoscopic 
time scale r^. 
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